
function [] = Chorin_Projection(dt)
global Ux Uy  h  g  p  q r dx dy L Nx Ny D E nt;


dt2= 1/dt;
beta=0.5;

dx1=1/dx;
dx2 = 1/(dx*dx);

dy1= 1/dy;
dy2 = 1/(dy*dy);

N2 = Nx*Ny;

%---- Construct the spatial Discretization Matrix------%
L=zeros(2*Nx*Ny,2*Nx*Ny);
A = generate_hyperbolic( Nx, Ny, Ux*0.5*dx1, Uy*0.5*dy1);
L(1:Nx*Ny,1:Nx*Ny)=A;
L(Nx*Ny+1:2*Nx*Ny,Nx*Ny+1:2*Nx*Ny)=A;
A = generate_elliptic(Nx,Ny,h*dx2,h*dy2);
L(1:Nx*Ny,Nx*Ny+1:2*Nx*Ny)=A;

for k=1:Nx*Ny
    L(Nx*Ny+k,k)=g;
end

A = generate_elliptic(Nx,Ny,h*h/3*dx2,h*h/3*dy2);
D=A;
C=zeros(2*Nx*Ny,Nx*Ny);
C(1:Nx*Ny,1:Nx*Ny)=A;

E = generate_elliptic(Nx,Ny,-2/15*h*h*h*dx2,-2/15*h*h*h*dy2);

for k=1:Nx*Ny
    E(k,k)=E(k,k)+h/3;
end
%---- Construct the spatial Discretization Matrix------%

I =eye(2*Nx*Ny,2*Nx*Ny);
r= zeros(Nx*Ny,1);

sol(1:Nx*Ny,1) =p;
sol(Nx*Ny+1:2*Nx*Ny,1) =q;

for n=2:nt+1;
    
    rhs = (dt2*I - (1-beta)*L)*sol;
    sol_new=(dt2*I+beta*L)\rhs;
    
    r=R_Solve(sol_new(N2 +1:2*N2 ));
    sol = sol_new + dt*C*r;
    
    p = sol(1:N2 ,1);
    q= sol(N2 +1:2*N2 ,1);    
    
%     if rem(n,100)==0
%         refreshdata(h1,'caller') % Evaluate p in the function workspace
%         drawnow
%     end
    
end
display('Completed Successfully');